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We implement the level set method for numerical simulation of the motion of a suspended particle 
convected by the fluid flow in a microchannel. The method automatically cope with the interactions 
between the particle and the channel walls. We apply the method in a study of particles moving in 
a channel with obstacles of different shapes. The generality of the method also makes it applicable 
for simulations of motion of particles under influence of external forces. 



I. INTRODUCTION 

In recent years numeral lab-on-a-chip systems have 
been developed to analyze biological samples. Many of 
these systems rely on handling of particles and cells com- 
parable in size to the dimensions of the channels contain- 
ing them. Examples of such microsystems are bumper- 
arrays or DEP-systems P, Q, 0, |j| 

It is a major challenge in theoretical microfluidics to 
study the dynamics of particles of finite size when they 
are convected by a fluid flow. Especially problematic is 
the forces appearing during collisions of the particles with 
the walls of the channel. 

The level set method [j| is well suited to cope with 
these problems. By introducing a hypersurface (/>(r,i), 
the particle interface is represented as the zero level set 
4>(r,t) = 0. The major advantage of the method is that 
this zero level set can be calculated implicitly instead of 
explicit tracking of the points on the interface. 

The manuscript is organized as follows: In Sec. HT1 we 
state the equations governing the dynamics of the system 
and in Sec. lIIII we derive the level set formulation for the 
tracked interface. The implementation of the method in 
the numerical simulation tool Femlab is described in 
Sec. lIVI and we present results of a test study in Sec. lVTI 
Finally, we evaluate the method in Sec. IVIII and give 
suggestions to future areas of usage. 



common boundary between fl\ and O2 is the interface T 
which we want to evolve. 

The rate of change of the momentum of the fluid is 
given by J n P%^ dr involving the substantial time deriva- 
tive of u. The change in momentum arises from the forces 
acting on the volume of fluid. In a microfluidic system 
we can neglect gravity and the only force F CT acting on a 
volume of fluid f2 stems from the stresses cr exerted by 
the surrounding liquid on the surface 9fi, 



F a = / cr ■ da, 



(2) 



where cr is the stress tensor modelled by 

&ij = -pS i: j + 7/ (d-jUi + diiij) . (3) 
Newton's second law therefore takes the form 



p— — dr = / er • da. 



(4) 



The right hand side of this equation can be split up in 
three integrals; two parts for each of the boundaries of 
the two subdomains and one along the common interface 



II, GOVERNING EQUATIONS 

We consider microfluidic systems. Hence the charac- 
teristic length scales of channels are of the order of 10 /im 
which is well beyond the intermolecular distances char- 
acteristic of the fluids involved. Thus the continuum hy- 
pothesis applies. Moreover, in these systems the flow 
velocities are much smaller than the propagation of pres- 
sure (the speed of sound) . We can therefore consider the 
fluids to be incompressible and the continuity condition 







(1) 



holds true for the velocity field u of the fluid. 

Consider a domain ft consisting of two subdomains f2i 
and Q.2 with surfaces dfli and dfl2, respectively. The 



p dr = / er-da+ / cr ■ da + / [cr • da 



Dt 



Oil, 



= V-erdr+ / V • er dr + / 7K da, 



(5) 



In the second equality we have used Gauss' theorem as 
well as the Young-Laplace law relating the pressure drop 
[cr ■ da] across the interface T to the surface tension 7 and 
average curvature k. 

To facilitate numerical computation it is desirable to 
rewrite the last integral in Eq. |H as a volume integral 
like the rest of the terms. This can be achieved by intro- 
ducing a level set function ^(r, t) as we will show in the 
following. 
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III. THE LEVEL SET METHOD 

Following Ref. Q we introduce a level set function 
(r, t) with the properties 



(r,i)>0, reOi, 
(r,t)=0, rer, 
(r,i)<0, re!! 2 . 



(6) 



This function uniquely defines the interface as T(t) = 
{r\(f>(r,t) = 0} and permits us to distinguish each sub- 
domain by the sign of <fi. We also introduce a transverse 
level set function tp(r, t) such that 



(7) 



We show in Appendix El that it is possible to construct 
such level set functions. In the following we consider a 
two dimensional system, but the method is applicable 
in higher dimensions also. We can construct a global 
orientation-preserving diffeomorphism that maps Q i— ► CI' 
through the variable change 



x' = tp(x,y) 
V = <P{x,y). 



(8a) 
(8b) 



We denote partial derivatives with indices, e.g., ip x = 
d x ip. The change of variables Eqs. <(BJ is area preserving 
because the Jacobian is non-zero, 



d{x,y) 



-^)-(^,^) = |V0||V^|^O ) (9) 



where we assume that ip is constructed such that Vt/> is 
parallel to the tangent direction and therefore — V^||V^. 

Furthermore we introduce a parameterization 
(x(s),y(s)) of r, where s is an arc-length variable. 
Using this parameterization an infinitesimal change in 
x' along T is given by 



dx% =0 = | VV'I ds, 



(10) 



where we have utilized the above assumption that the 
gradient of -0 is parallel to the tangent direction. With 
the above definitions we can rewrite the surface integral 
in Eq. © as 



7«da= / 7Knds 
I J<p=o 



=o 7K |V0||W| 

sr/ n 1 
n , \V(j>\ \Vlp\ 



dx' 



(11) 



dx'dy' 



where we have used that the normal n to the interface 
can be written as V0/|V0|. Using Eq. 10 for changing 
variables, Eq. ill 111 becomes 

/7Kda= / -fK5{<j))V^dxAy. (12) 



Inserting Eq. ltT2*|l into Eq. © yields 

/ p^dr = ( [V • a + 7K5(4>) V<j>] dr. (13) 



This must hold true for any volume ft. Hence 



p [d t u + (u • V)u] 



(14) 



which is the level set formulation of the Navier-Stokes 
equation. 

In order to have the system completely described by 
dynamical equations we finally need an equation describ- 
ing the evolution of the zero level set. We only need to 
consider the movement of the zero level set because this 
is the only part of the level set function with a physical 
interpretation. Evolving the equation </>(r, t) = in time 
defines the movement of the front. Differentiating with 
respect to time yields gr0(r, t) — which is written as 



where V 



dr 

dt 



a^(r,t)+V-V0(r,t)=O, (15) 
is the velocity of the zero level set. 



Requiring the velocity field to be continuous leads to 
V = u, and the evolution equation for <f> becomes 



0. 



(16) 



IV. FEMLAB IMPLEMENTATION 



One of the great advantages of the level set formula- 
tion is that it does not track the interface explicitly but 
rather capture it implicitly. Thereby we avoid to intro- 
duce explicit forces from the walls during collisions as 
they enter implicitly through the stress tensor er and the 
no-slip boundary condition on the velocity field u. Fur- 
thermore, several numerical tools are available for solving 
the dynamical system. In this section we describe how 
we have implemented the level set method in the finite 
element software package FEMLAB Q. We have used the 
Femlab scripting language trough a Matlab interface 
in the general PDE mode. Here the PDEs are given by 



d a - 



dU 

dT 



v r = f 



in Q, 



(17a) 



in terms of the variable vector U, the current tensor T 
and the generalized source terfield F. The boundary con- 
ditions take the form 



OU, 



= R„ 



on d£l 
on 90, 



(17b) 
(17c) 



where the index I is the variable counter, m is the con- 
straint number (the number of boundaries) and j is the 
number space dimension number. The Lagrange multi- 
pliers p m are chosen by Femlab in order to fulfill the 
constraints, while the scalars Ft, Gi and R m are given by 
the physics of the problem. 
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A. Navier Stokes equation in Femlab 

Introducing the characteristic length scale Lq, velocity 
scale U , density po, viscosity t] and surface tension 70 
we can express the physical quantities as a dimension- 
less number times the characteristic scale. Denoting the 
nondimensional quantities by a tilde we simply have 



Lqt, 
VoV, 



u 

7 



U u, 
7o7- 



P = PoP-, 



(18) 



Similarly we can define the characteristic pressure and 
timescale as relations between the chosen characteristic 
parameters 



P : 



V0U0 _ 



Lo- 
ft = —t, 
U 



(19) 



Substituting Eqs. iflSl) and ltl9|l into the Navier-Stokes 
equation ifHj) yields 



Rep dfU + (u • V)u 



Ca 



jK8((f>)V<f>. (20) 



Here the Reynolds number Re = p UoL /r]o is the ratio 
between inertial forces and viscous forces and the Capil- 
lary number Ca — rjoUo/jo is the ratio between viscous 
forces and the surface tension forces. 

Rearranging the terms in Eq. l|20j) we find 

1 



i?e/5<9 t ~u — V • cr = -—jft5((f))'V( l 



i?ep(u-V)u, (21) 
which is seen to be on the Femlab general form if 



d a = Rep, 

r = -ar. 



F = -Rep(u ■ V)u + —jkSU) V6, 
La 



U u =u. 



(22a) 
(22b) 

(22c) 
(22d) 



The density p, viscosity 77 and the curvature of the front k 
are defined as auxiliary functions of the level set function 
4>. In a system with two immiscible incompressible fluids 
(or a particle in a fluid) the density and viscosity are 
constant on each side of the interface. We can therefore 
define the dimensionless density and viscosity as 



and 



p=l + H(cf>) [ * - 1 

\P2 



1=1 + HW\£-1 
where H{cf>) is a Heaviside function defined as 

VW )0, 0£fl 2 . 



(23) 



(24) 



(25) 



Setting po = pi ensures that the density of the fluid is 
pi and pi in fix and 2 , respectively. Similarly setting 
Tjo = V2 makes the viscosity of the fluid r\\ and 772 in fii 
and VI2, respectively. 

The curvature of the zero level set is given by 



= V • n = V 



(26) 



where n = V0/|V^>| is a unit normal vector to the in- 
terface nn. 

When solving the system numerically the abrupt 
change in density and viscosity across the interface causes 
numerical instabilities to occur. In order to avoid this we 
substitute H(4>), d((j>) and sign(0) with the smeared out 
versions H e (<fi), 6 e (4>) and sign £ (</>) defined as 



5 e ((f>) = H' e ((p) = ^ - — tanh z 



1 1 

Ye ~ Yt 



(27a) 
(27b) 

sign e (0) = tanh . (27c) 

This implies that the interface has a finite thickness r e 
approximately given by 



2c 



(28) 



B. The continuity equation in Femlab 

The dimensionless form of the continuity equation is 
= V ■ u, (29) 

which is entered into Femlab by choosing F = V ■ u, 
r = 0, d a = and U p = p. 

C. The level set equation in Femlab 

The nondimensionalized form of the convection equa- 
tion for the zero level set is 



4>- t + u • Vcb = 0, 
which can be rearranged to 

61 = — u • W<h 



(30) 



(31) 



and implemented in Femlab by setting F = u • V0, 
r = 0, d a = 1 and Ua, — 4>. 



TABLE I: The parameter values used in the simulation of the 

test case. 



Reynolds number 


Re = 


1 x 10" 3 


Capillary number 


Ca = 


1 x 10 6 


Density 




1 x 10 3 kg m -3 


Viscosity 


'/U 


1 x 10" 1 Pas 


Obstacle size 


l = 


6 x 10" 6 m 


Particle radius 


r p 


3 x 10" 6 m 


Pressure drop 


Ap = 


1.2 x 10" 3 Pa 


Time step 


At = 


5 x 10~ 2 s 


Mesh element size 


hmesh = 


1.1 x 10~ 6 m 


Thickness parameter 


e = 


0.5 X hmesh 



D. Reinitialization of the level set function 

It is necessary to maintain a uniform thickness of the 
interface throughout the calculations. This requires that 
the gradient of the level set function is constant within 
a region around the interface < e. This is not au- 
tomatically fulfilled. The time evolution of any level set 
4>(r,t) = C is given by the level set Eq. itTfi)) . This means 
that the height of the level set function will remain con- 
stant, but it does not ensure that the gradient does not 
change. Thus in order to keep a fixed interface thick- 
ness we need to reinitialize the level set function without 
changing the zero level set. 

In principle we can use any function that fulfills 
Eq. ©, since only the zero level set has a physical in- 
terpretation. But requiring the interface thickness to be 
fixed constrains the gradient of (f) to be fixed in a region 
around the interface. A choice of <fi(r, t) that fulfills these 
requirements is the signed distance function, where the 
distance is the shortest distance d(r) from a point to the 
interface 

d(r) =±min(|r-r r |), (32) 

IT being the points on the interface. The plus sign applies 
if r G fix and the minus sign if r S O2 . The length of the 
gradient for this particular choice of level set function is 

|V0| - 1. (33) 

We have implemented two different reinitilization pro- 
cedures. One simple reinitialization procedure where we 
recalculate the level set function at every time step and 
one using the reinitialization equation suggested by Suss- 
mann, Smereka and Osher 0] 

d T ip(r,r) = signal - |V^(r,r)|), (34) 

with the initial condition t{i(r, 0) = 4> and r being a pseu- 
dotime. The steady state solution to this equation is the 
reinitialized level set function. Because numerical oscil- 
lations can occur if the sign of <\> changes abruptly at 
the interface it is necessary to use the smeared out sign 
function given in Eq. i|27c|) . 
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FIG. 1: For the test study we use the geometry and mesh 
shown in the figure. The general shape of the obstacle is as 
shown in the lower inset on the right. The radius a of the 
rounded corner was changed from one simulation to the next. 
The aspect size of the obstacle is The height of the channel 
is H = (20/3)/ and the width of the channel is W = (13/3)/. 
The upper inset on the right shows the general idea of the 
test study: The particles start in the initial position xo and 
the final position xansd is recorded. 

The reinitialization equation is already on a form suit- 
able for implementation in Femlab. Simply letting F 
equal the right hand side of the equation and setting 
d a = 1 and r = with = ip does the trick. 

To avoid mass loss during the reinitialisation procedure 
we have put a constraint on the solution: the volume of 
the particle must be constant at all time. This is done 
in Femlab via the field f em. equ. constr where we con- 
strain the difference between the integrals of the smeared 
out Heaviside function H e (tp) at time r and the smeared 
out Heaviside function -ff e (</>) at time t = to be zero. 
The integrals are computed by using the integration cou- 
pling variables in Femlab. 

V. MODEL SYSTEM AND SETUP 

To test the implementation of the level set method in 
Femlab we have done a test study of a particle (a drop 
of high viscosity and surface tension) which is passively 
convected in a two dimensional fluid flow. The viscosity 
i]2 of the particle was 100 times larger than the viscosity 
771 of the fluid. The density p\ of the fluid was equal 
to the density P2 of the particle. The complete list of 
parameters is given in Table Q] 

The physical domain is an infinitely wide and infinitely 
long channel with an obstacle in the center as shown in 
Fig-El The boundary conditions on the fluid are no-stress 
on the sides of the computational domain and no-slip at 
the obstacle. The fluid velocity field is periodic from 
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0.2 0.4 0.6 0.8 1 

2a /I 

FIG. 2: For particles passing obstacles of different shapes nor- 
malized difference 2/S.x/W in horizontal position from start 
to finish is plotted versus starting position 2a//. The missing 
data points for the simulations with the initial positions of 
the particles nearest to the center of the channel is due to the 
particles getting stuck at the obstacle and hence not reaching 
the final position. 



top to bottom of the domain and is driven by a pressure 
difference Ap. 

We ran a series of simulations with the shape of the 
obstacle changing from circular to quadratic by changing 
the radius of the rounded obstacle corner a. Each sim- 
ulation consisted of a series of runs with different initial 
horizontal position xo of the particles and the initial ver- 
tical position of the particles was yo = H — l from the top 
of the channel. When the center of a convected particle 
is I from the bottom of the channel the final horizontal 
position Xfl na i is detected (Fig.QJ. 

We represent the particle by the negative part of a level 
set function and the surrounding fluid is identified by the 
positive part of the level set function. The initial level 
set function is given by 



<fr(x, y, t = 0) = ^(x - x ) 2 + {y- yo) 2 - r p , (35) 

where (xo,yo) is the initial position of the particle and 
r p is the radius of the particle. Using these parame- 
ters we solve the problem by first evolving the dynamical 
equations in a small time step At and then reinitialize 
the level set function using the reinitialization procedures 
described above. With the reinitialized level set function 
as initial condition for (j> we evolve the dynamical system 
one more time step. This sequence is continued until the 
particle has moved all the way through the system. 




x 



final 



FIG. 3: The paths of particles passing obstacles of different 
shapes when the starting point is 2xo/W = 0.308 right of the 
centerline of the channel. 




FIG. 4: The path of the particle started at 2x /W = 0.015 
when the radius of the rounded obstacle corner is a — 1/2. 
The particle (black dot) is shown when it 'interacts' with the 
obstacle. The small gap between the particle and the obstacle 
wall is caused by the smearing of the particle interface. 



VI. RESULTS 

We carried out simulations for four different initial 
positions of the particle. The initial horizontal posi- 
tions 2x /W were 0.015, 0.077, 0.308 and 0.539, re- 
spectively. For each of these initial positions we used 
five different radii of the rounded corner of the obstacle: 
2a/l = i/10, with i = 1,3, 5, 7, 10. 

For each combination of initial position and obstacle 
shape we solved the system and obtained the particle 
paths. Examples are shown in Figs. and 31 It is seen 
that the paths of particles with the same initial position 
changes as function of the shape of the obstacle (Fig. El . 
In Fig. we have plotted the difference in the horizontal 
position Ax from start to finish. 

The difference in horizontal position is almost zero for 
the particles started in at the greatest distance from the 
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center of the channel, independent of the shape of the ob- 
stacle. As the initial position gets closer to the center of 
the channel the difference in horizontal position becomes 
larger and the round obstacles tend not to drag as much 
in the particles as the square obstacles yielding a larger 
difference in the horizontal position. 

Fig. Q] shows that our implementation of the level set 
method is capable of coping with the interaction forces 
between the stable obstacles and the moving particles 
automatically. 



VII. DISCUSSION AND CONCLUSIONS 

We have shown that the level set method is easily im- 
plementable in Femlab and that it is a suitable method 
for coping with the interaction forces between particles 
and hard walls automatically. Particles can be modelled 
as very viscous liquid drops and the shape preservation 
can be taken care of trough an appropriate reinitializa- 
tion procedure. 

We have used a simple shape preserving reinitializa- 
tion method. Further work is needed in order to con- 
vert particles of an arbitrary fixed shape. One promising 
reinitialisation scheme is t he p article level set method 
suggested by Enright et al. [10] . 

The level set method might prove useful when simu- 
lating microfluidic systems for particle handling. In this 
paper we have only considered the forces exerted on the 
particles by the convecting fluid and thereby indirectly 
the forces from the solid walls. However also other forces 
such as DEP forces or magnetic forces could be taken into 
account making the method applicable for simulations of 
many lab-on-a-chip systems fabricated today. 



APPENDIX A 

We demonstrate how to construct the transverse level 
set function ip with the required properties. We start by 
defining a coordinate transformation by 

^- (x(s, r),y(s, r)) = V(p(x(s, r), y(s, r)) , (Ala) 
where 

(x(s,0),y(s,0^ = (x(s),y(s)y (Alb) 

Because of the 5 function in Eq. itTHl ip only needs to 
fulfill the requirements in a small region |t| < e around 
T. In this small region we can define ip as 



(A2) 



where ipo ( s ) is a smooth increasing function if and only 
if the mapping of (x, y) to (s, r) is one-to-one. Using the 



change of variables theorem we have to show that 
d(x,y) 



(A3) 



d(s,r) 

Taylor expanding Eq. IIAla]) around r = yields 

(x T ,y T ) = V^(x(s),y( S ))+0(r). (A4) 



Differentiation of Eq. ijAlall with respect to s and inte- 
gration with respect to r yields 



^ T ^V0(x( S ,T'),y( S ,r'))dr' (A5) 
From which follows 

(x s (s,T),y(s,Tfj - (x s (s,0),y s (s,0)) = 

^ T ^V0(.T( S ,r'),2/( S ,r'))dr', (A6) 

and thus 
(x s (s,T),y(s,T)\ 

= (x a (s),y B (s)) + J —Vct>(x(s,T'),y(s,r'y)dT' 
= T(s) + 0(t). (A7) 

Here T is a unit tangent vector to the interface. We can 
now calculate the determinant i|A3j) 



d(x,y) 



= (x T ,yr) ■ (-y s ,x T ) 

= Vcf>(x s> y s )-f 

= |V0||T|+O(r) 

= |V^ =0 + O(r)^0. 



(A8) 



This means that ip is well defined in a small region around 
r. Now all we need to prove is that and VV> are 
orthogonal and that \Vtp\ ^ 0. The orthogonality can be 
proved by differentiating -0 with respect to r, 



— tp(x(s,T),y(s,rjj = ip x x T + ^ v y T 



dr 



(A9) 



which means that (ft and ip are orthogonal if and only if 
|Vi/'| ^ 0. This follows immediately from differentiating 
ip with respect to s, 



-il)(x(s,T),y(s,Tfj = ijj x x s +ip y y s 



Vip ■ (x s ,y s ) 
Vip ■ T 

|V^|=^(s)>0, 



(A10) 
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because ipo(s) was chosen to be an increasing function. of the Navier-Stokes equation for a two liquid flow of 
Thereby we have established the level set formulation incompressible fluids. 
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